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We investigate the effect of single and multiple impurities on the Zeeman-localized, spin polarized 
bound states in dilute magnetic semiconductor hybrid system. Such bound states appear whenever 
a dilute magnetic semiconductor showing giant Zeeman effect is exposed to an external magnetic 
field showing nanoscale inhomogeneity. We consider the specific example of a superconductor-dilute 
magnetic semiconductor hybrid, calculate the energy spectrum and the wave functions of the bound 
states in the presence of a single impurity, and monitor the evolution of the bound state as a 
function of the impurity strength and impurity location with respect to the center of the Zeeman 
trapping potential. Our results have important experimental implications as they predict robust 
spin textures even for than than ideal samples. We find that for all realistic impurity strengths the 
Zeeman bound state survives the presence of the impurity. We also investigate the effect of a large 
number of impurities and perform ensemble averages with respect to the impurity locations. We find 
that the spin polarized Zeeman bound states are very robust, and they remain bound to the external 
field inhomogeneity throughout the experimentally relevant region of impurity concentration and 
scattering strength. 



I. INTRODUCTION 

Due to its potential use in spintronics and quantum 
information applications, the controlled manipulation of 
individual spins has been a subject of great interest in 
recent years. In the search for appropriate materials for 
device fabrication, one possible route is the use of III-V 
or II- VI diluted magnetic semiconductors (DMS). 

In II- VI DMSs, the exchange interaction between band 
charge carriers and the localized magnetic spins of the 
magnetic dopant (such as Mn) gives rise to a giant split- 
ting between band states with different spin components. 
This giant Zeeman spin splitting has been extensively in- 
vestigated and characterized^ Given the linear depen- 
dence of this splitting on the external applied magnetic 
field for fields lower than roughly IT, it can be described 
in terms of a very large effective g-factor for the band 
carriers. For example, Dietl et al. 3 reported an electron 
effective g-factor of about 500 at sub-Kelvin tempera- 
tures in a CdMnSe sample, which implies a value of about 
2000 for the effective ^-factor of a hole in this material. 
In previous works^ 7 -^ we proposed to take advantage 
of this large effective ^-factor by combining it with a spa- 
tially inhomogeneous applied magnetic field, to create a 
spatially varying Zeeman potential that acts as a confin- 
ing potential for only one spin orientation. In suitable 
conditions, this leads to a single charge carrier with a 



well defined spin being trapped in the region of large lo- 
cal magnetic field. This spin-polarized charged carrier 
can then be manipulated through external control on the 
applied inhomogeneous magnetic field. 

The needed nanoscale inhomogeneous magnetic field 
can be generated in various ways. One possibility is to 
use nanoscale magnets placed above a DMS quantum well 
(QW) ^ 6 ' 7 ' 10 ' 11 Depending on the shape and orientation 
of the nanomagnet, different nonuniform fields are gen- 
erated, giving rise to various types of confined states^ 

Another possibility is the use of Abrikosov vortices 
that appear in type-II superconductors (SCs). Above the 
lower critical field B Cl vortices populate the superconduc- 
tor, forming a vortex lattice. Nano-engineering can be 
used to insure that the vortices nucleate in well defined 
positions^ The field of a single vortex is non-uniformly 
distributed around a core of radius r ~ £ (£ is the SC 
coherence length) and decays away from its maximum 
value at the vortex center over a length scale A (A is the 
penetration depth) . If a SC layer that hosts such vortices 
is placed above a DMS layer (QW), the inhomogeneous 
magnetic field of the SC vortices creates an inhomoge- 
neous magnetic field in the DMS layer. According to our 
previous calculations^ these fields are sufficiently large 
to result in the confinement of band carriers with a given 
spin orientation in the small region of the DMS QW that 
is located directly under a SC vortex core. Thus, spin 
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textures are formed and can be manipulated by moving 
the source of the magnetic field, i.e. the SC vortex. In 
other words, the SC vortices act as spin and charge tweez- 
ers and can be used for a wide array of applications, from 
investigation of the Hofstadter butterfly and making spin 
shuttles^ to generating anions of interest for topological 
quantum computing^ 

In our previous work we assumed that the DMS QW 
is perfectly clean and with Mn impurities distributed in 
a perfectly homogeneous fashion (for instance by digi- 
tal doping^), so that the Zeeman potential is a smooth 
function mirroring the applied magnetic field. For such 
systems we obtain typical binding energies of the spin- 
polarized charge carrier on the order of a few meV4^ As 
we show in this paper, such perfect homogeneity is not 
necessary. This is a very important result from experi- 
mental detection and application point of view, since it 
opens up the possibility of investigating Zeeman local- 
ization effects in samples grown during standard molec- 
ular beam epitaxy (MBE) conditions. Since the effec- 
tive (/-factor is proportional to the local concentration of 
Mn, inhomogeneities in their distribution will result in 
a spatially- varying g(r) which will modify the trapping 
potential. Even more detrimental are charged impurities 
which exist in any sample. Charge dopants are needed to 
introduce charge carriers in the II- VI DMS, because the 
Mn is isovalent with the valence-II element it replaces. 
The disorder due to these charged dopants can be min- 
imized in the usual way used for two-dimensional elec- 
tron systems, by doping at some distance away from the 
DMS QW and using gates to control the concentration of 
charge carriers in the DMS QW. However, small concen- 
trations of undesired charged dopants in the DMS QW 
cannot be avoided. 

In this article we investigate the role played by charged 
scattering centers on the confinement of spin polarized 
carriers in the DMS QW. We focus on repulsive centers, 
which are expected to be most detrimental to the bind- 
ing of the spin-polarized carrier in the Zeeman trap, and 
then briefly discuss the other types of scattering centers 
mentioned above. We show that due to the difference 
between the length scale of the confining potential and 
the repulsive scattering potentials, the bound state en- 
ergy is not significantly changed even in the presence of a 
large number of impurities. This result indicates that the 
binding of spin-polarized charge carriers in Zeeman traps 
created in a DMS QW is very robust against such defects, 
with immediate implications for spintronics applications 
based on this scenario. 

It is important to also emphasize that our conclusions 
apply to a much wider range of systems than the one 
of interest to us, because these conclusions are based on 
quite simple physics, as we show below. We argue that in 
any system (such as quantum dots, quantum wells, etc.) 
where localization occurs on a lenghtscale much larger 
than that of typical scattering potentials, the effect of 
such scattering centers is very small even if their potential 
is very large. 



This paper is organized as follows: In the next section 
we describe in more detail the system we study, our the- 
oretical approach and the approximations we have made. 
In order to build intuition and to help interpret later cal- 
culations, in section HTTl we analyze results for a simplified 
one dimensional (ID) model in the presence of a single 
impurity. This allows us both to gauge our computa- 
tional scheme against exact solutions, and to gain some 
intuition about the effect of the scattering centers. In 
section HVl we repeat the analysis for a single scattering 
center for the two-dimensional system of interest to us, 
and discuss the influence of the symmetries in our results. 
Finally, in section [V] we present results for a random dis- 
tribution of impurities and various impurity concentra- 
tions. We conclude in section IVT1 with a summary and 
discussion. 



II. 



THEORETICAL MODEL 



The hybrid structure of interest to us is sketched in 
Fig. Q] ( a ) > and consists of a SC layer in the vortex phase 
placed above a DMS QW. The two materials are sep- 
arated by a thin insulating layer. The superconductor 
can be Pb or Ni, for example - both these metals can be 
grown on top of semiconductors using MBE techniques. 16 
The QW consists of a weakly p-doped DMS. An example 
of such a DMS QW is a p-doped (Cd,Mn)Te well, doped 
with N by a modulation doping technique that avoids 
inhomogeneity effects and increases the mobility of the 
carriers^ The carriers of interest to us are confined in 
the 2D DMS QW. The giant Zeeman effect in the DMS 
is described by a very large g-factor. This combines with 
the inhomogeneous magnetic field generated by each vor- 
tex to give rise to an effective potential that binds a spin 
polarized hole under each vortex, in a clean sample. As 
stated, we are now interested in the effect of repulsive 
scattering centers on these spin-polarized bounds states. 
A typical potential well in the presence of a single, and 
of several charged scatterers is illustrated in Figs. Gib) 
and[ljc), respectively. 

In a parabolic approximation, the Hamiltonian for a 
charge carrier in the 2D DMS quantum well which in- 
teracts with the magnetic field of the vortex and with 
scattering centers is 



H = 



[p - qA{r)f 
2m 



^ffoff Mb*? • B(r) 



Vi(f), (1) 



where m and q are the effective mass and charge of the 
carrier, B{f) is the magnetic field generated by a sin- 
gle vortex, A(r) is its vector potential and Vj(f) are the 
repulsive potentials from the i scattering centers. For 
simplicity, we assume that the motion of the carrier in 
the QW is two dimensional. In the calculations shown 
here, we use m = 0.5m e and g e s =500, these being typical 
values for holes in II- VI semiconductors. Because of the 
large g c s, for a single vortex the effect of the A(x, y) term 
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FIG. 1: (color online) (a) Illustration of the magnetic field 
generated by a single vortex in a hybrid system made of a 
type II superconductor on top of a DMS quantum well which 
contains impurities, (b) Potential profile of the Zeeman trap 
plus a single repulsive impurity and (c) same, but for a ran- 
dom distribution of impurities (figures not to scale). 



is negligible compared to that of the Zeeman interaction 
and we ignore it from now on£ 

Experimental results show that the Zeeman splitting 
for holes in a DMS is anisotropic, depending on the di- 
rection of the magnetic field with respect to film plane. 17 
Using a Luttinger Hamiltonian^ one can also show that 
in a 2D QW, the g e ff of heavy holes is highly anisotropic, 
with an in-plane component much smaller than that per- 
pendicular to the film plane. 7 With these results in mind, 
we can further simplify our problem by considering only 
the effect of the z component of the magnetic field, which 
is most strongly coupled to the charge carrier. This ap- 
proximation allows us to decouple the Hamiltonian in the 
spin-up and spin-down sector, and consider them sepa- 
rately. We will only focus on the spin-component for 
which this Zeeman potential acts as a trap. 

As discussed extensively in Ref. H, the size of the 
confining region defined by the inhomogeneous magnetic 
field depends on the superconductor's parameters and the 
distance between the SC and DMS layers. It is typically 
several tens of nm in size. On the other hand, the im- 
purity potential of the scattering centers is considerable 
only in a much smaller range of a few A in the immediate 
vicinity of the impurity. It is the effect of this strong re- 
pulsive potential on the bound state that is of interest to 
us, since the long-range Coulomb repulsion is too weak 
to be relevant (although its effects can also be studied 
with our method, if so desired). Given the large differ- 
ence between the characteristic length scales, we model 
the scattering potential as a delta function. 

After all these approximations, the Schrodinger equa- 
tion for the trapped spin component is: 
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where a characterizes the strength of the impurity poten- 



tial and Ri is the location of the i th impurity. 

We solve for the eigenvalues in the usual way, by ex- 
panding the wavefunction in a complete basis set and 
finding the appropriate coefficients from a matrix equa- 
tion. Our basis functions are B-spline polynomials on a 
non-uniform knot sequence adjusted optimally for each 
specific arrangement of charged impurities. This method 
has been widely used in atomic physics^ and is ideally 
suited for problems such as this. Details of the method 
and its advantages are discussed in the appendix [A] 



III. ONE DIMENSIONAL MODEL 

We first study a one dimensional model with a single 
delta function impurity inside a square potential well. 
This provides an intuitive understanding of the physics 
we want to explore, and also a numerical test of our B- 
spline scheme since it permits the comparison between 
the numerical results and the exact solution. 

We consider a ID square well potential 



V (x) = 



-do, if - L/2 <x< L/2 
0, otherwise 



(3) 



where L is the lateral size of the well, 
single impurity potential 

V\{x) — a8{x — xq) 

where xq is the impurity position. 
Schrodinger equation is given by 



To it, we add a 



(4) 

In this case, the 



2m dx 2 



ip(x) + V (x) + Vi(x) = eip(x) 



(5) 



To understand the effects of the impurity on the bound 
states for different strengths of the impurity potential, 
we calculated the eigenvalues of this Hamiltonian. The 
energy of the lowest two bound eigenstates is shown in 
Fig. [5] vs. the strength of the impurity potential ex- 
pressed as a dimensionless quantity a = a/(aoL). The 
parameters are L = lOOnm and ao — 3.5ft 2 /(mL 2 ), and 
the impurity is located in the center of the well, xq = 0. 
Fig. [2] reveals that the ground-state energy first increases 
with a and then saturates to the value of the first ex- 
cited state. On the other hand, the energy of the first 
excited state is independent of a. The reason for this be- 
havior is straightforward to understand. Since the first 
excited state has a node at the location of the impu- 
rity, its energy and wavefunction is naturally unaffacted 
by its presence. On the other hand, the original ground 
state wavefunction has a maximum at the impurity loca- 
tion, so its energy increases with a. However, for large 
enough a the ground-state wavefunction also develops a 
dip precisely at the location of the impurity and becomes 
insensitive to its strength. This is demonstrated by the 
insets to FigJ31 which show the ground-state wavefunc- 
tion at a = and a = 100. In essence, when a > 1, 
the impurity potential splits the original well into two 
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0.0001 0.01 1 100 



a 

FIG. 2: The ground state and first excited state energy versus 
a, a dimensionless quantity related to strength of the impurity 
potential. The impurity is located at the center of the well 
x — 0. The left inset shows the ground-state wavefunction for 
a = 0, and the right inset shows it for a — 100. 



isolated wells which are degenerate, given the symmetry 
of the problem. 

When the impurity is placed off-center, as in Fig. O the 
system loses this symmetry and the 8 function separates 
the original well into two distinct wells of different sizes. 
Since the impurity is now placed in a location where both 
the original ground state and the first excited state wave- 
function were finite, both energies now increase with a 
when a < 1. At higher values of a both the ground state 
and excited state energies saturate. As shown in the in- 
sets, the wavefunctions develop dips at the location of 
the impurity and become progressively more localized on 
one or the other smaller wells. The new ground state 
wavefunction becomes confined to the larger new well. 

While these results are very easy to understand, they 
do provide key insight into the effects of the impurity on 
the bound states. First, it is clear that symmetry plays 
a very important role. While this influence is maximized 
by the one dimensional character of the model, we expect 
analogous features in the two-dimensional case. 



- 




0.01 1 100 
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FIG. 3: The ground state and first excited state energy versus 
the strength of the delta function impurity for xo = L/4. a — 
a/aoL. The left inset shows the ground-state wavefunction 
for q = 0, and the right inset shows it for a = 100. 



More importantly, we already see that even if the im- 
purity potential is large, its overall effect on the bound- 
state energy and wavefunction is limited. Because it is 
extended over a much larger area, the bound wavefunc- 
tion simply avoids the impurity by developing a dip (or 
a zero) in its vicinity, at relatively low energetic cost. 
In ID the impurity effectively "cuts" the wavefunction. 
This will not happen in d = 2 or higher dimensions. How- 
ever, as we demonstrate below, the limited effect on the 
binding energy is a general feature observed in higher 
dimensions. 

The results shown so far also illustrate why B-splincs 
are an ideal basis for such problems. If one tried to 
use a more traditional basis of orthonormal polynomials 
such as harmonic oscillator eigenstates, one would have 
to mix in very many basis states in order to be able to 
accurately describe smooth wave-function variations on 
a short length scale, near the impurity, as shown in the 
right inset of Fig. [51 Describing regions where the wave- 
function essentially vanishes, as shown in the right inset 
of Fig. [31 would be even more difficult. On the other 
hand, use of B-splines allows one to just slightly increase 
the number of basis states by sampling the vicinity of the 
impurity on a smaller mesh. As a result, the calculation 
with one or more impurities is very comparable, in terms 
of numerical computational costs, with the one in the ab- 
sence of impurities. A more detailed discussion of these 
issues is given in Appendix A. 

IV. TWO DIMENSIONAL ZEEMAN 
POTENTIAL WITH A SINGLE IMPURITY 

We now can consider the two dimensional system in 
the presence of a single charged impurity. As we are 
not interested in the details of the Zeeman potential, we 
take the magnetic field generated by the vortex to have a 
Gaussian profile. This approximation is reasonable close 
to the vortex core, where the magnetic field decays ex- 
ponentially on a length scale set by £j2i2£ We therefore 
solve Eq. [2] where the magnetic field profile is taken to 
be: 

B z (r) = B exp (-(r - n,) 2 /? 2 ) , (6) 

where ro is the location of the vortex core, and Bo = 
0.206T and £ = 50 nm are the strength and the range of 
the magnetic field (these are typical values for a dirty Pb 
film in a type II regime). We solve the equation numer- 
ically for a system of finite size 200nm x 200nm, with 
periodic boundary condition and a grid of 51 x 51 knots. 
This large lateral size is chosen so as to eliminate finite 
size effects in our bound state energies. For more details 
about the implementation of this numerical method and 
the periodic boundary conditions, see the appendix. 

We follow a procedure similar to the one in the pre- 
vious section and begin with a single impurity located 
at the center of the Gaussian and vary its potential 
strength. Again, we define a dimensionless quantity 
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FIG. 4: (color online) The two lowest state energies versus 
the impurity strength. The left(right) inset is the ground- 
state wavefunction for a — ( a = 1). 

a = 2a/(g c fffiBBoirt; 2 ) related to the strength of the im- 
purity potential. When the impurity is located at ro, as 
shown in Fig. 0] the ground state energy increases with 
the a and saturates at a finite value for a > 0.01. In the 
insets of Fig. [H we can see that the ground-state wave- 
function maintains its s-type symmetry as we increase a 
but develops a dip at its center. Similarly to the one di- 
mensional case, the first excited state, which is two-fold 
degenerate, is not affected by the impurity potential, as 
can be seen in Fig. [5] (a) and (c). This is because its 
eigenf unctions have p-type symmetry, with a node at the 
origin where the impurity is located. These results mir- 
ror those of the ID case but are less dramatic: the en- 
ergy separation AE between the two lowest eigenstates 
changes by less than 30% for arbitrarily strong repulsive 
impurity potentials even when the charged impurity is 
placed where it should do the most damage. 

Next, we investigate the dependence of the bound en- 
ergies and wave-functions on the location of the impu- 
rity, for a fixed strength a. The results are shown in 




FIG. 5: (color online) Wave- functions of the two excited states 
with p-type symmetry when the impurity is located at the 
center of the vortex, in (a) and (c); and for an off-center 
impurity, in (b) and (b). 




distance film) 



FIG. 6: (color online) The lowest eigenenergies versus the dis- 
tance of the impurity from the center of the vortex. The solid 
line is the ground state energy, the dashed line is the energy 
of the excited state whose wavefunction is shown in Fig. [5] 
while the other excited state's wavefunction is shown in Fig. 
0(d). The left (right) inset is the ground-state wavefunction 
for n — (ri — 5 nm). 

Fig|6] The ground state energy decreases monotonically 
as the impurity moves away from the vortex center be- 
cause it costs less energy to create a dip at its location 
(the ground-state wavefunction is shown in the inset). 
Once the impurity is farther away than the character- 
istic length scale of the bound state, the ground-state 
energy saturates to its unperturbed value. Due to the 
symmetry of the problem, the degeneracy of the first ex- 
cites state is lifted. The state which has its nodal line 
where the impurity is located is unaffacted by it. The 
other eigenstate has a finite wavefunction at the impu- 
rity location, and its energy is increased. However, if the 
impurity is placed far enough, its effect vanishes and the 
two excited eigenstates become degenerate again. The 
excited wavefunctions are shown in Figs[5jb) and (d). 

Finally, we analyze how the energies of the lowest 
eigenstates change when the impurity is off-center at 
a fixed distance = 5 nm and we vary the impurity 
strength. The results are shown in Fig. [7j As expected, 
the ground-state energy increases but very little, since 
for large a a dip appears at the location of the impurity 
and its effect saturates, irrespective of its strength. The 
excited state with a nodal line at the impurity location 
continues to be unaffected by it, while the other excited 
state behaves similarly to the ground-state: its energy 
increases with a, but only by a finite amount until its 
wavefunction acquires a zero where the impurity is. 

These results clearly demonstrate that a single impu- 
rity is unable to unbind the Zeeman trapped charge car- 
rier, irrespective of how large its repulsive potential is. 
All it can do is to simply "poke a hole" in the wavefunc- 
tion, and this costs a relatively small energy. Even if the 
impurity is placed at the center of the Zeeman potential 
well, the energetic cost is little. This suggests that one 
could have many impurities in the vortex area without 
affecting the bound state significantly. Our expectations 




FIG. 7: The ground, first and second excited state energy 
versus the impurity strength. The left inset is the wavefunc- 
tion for a = 0, and the right inset is the wavefunction for 
a = 0.01. 

are indeed verified in the next section. 



V. RANDOM DISTRIBUTION OF IMPURITIES 

We now address a more realistic situation of a random 
distribution of multiple impurities inside the QW. In this 
case, as detailed in the appendix, the use of B-splines is 
particularly useful to minimize computational costs. 

For a large number of impurities, the energy and wave 
function will depend on the distribution of impurities. 
We follow the standard procedure and average our re- 
sults over various impurity configurations. We performed 
averages over 1000 configurations at different impurity 
concentrations. The digitally doped Cdi-^Mn^Te has 
a Zinc-Blende structure with a lattice constant of about 
6.5 A. There are 4 Cd in an unit cell, and we take x ~ 1.6. 
Assuming a single layer of Cdi-^Mn^Te gives a Mn areal 
concentration of 0.15 mn -2 . As already discussed, the 
number of charged dopants can be minimized by the use 
of co-doping techniques. While is impossible to eliminate 
all these scattering centers, we expect that their concen- 
tration is much smaller than the Mn concentration. 

Fig. [5]shows the central result of this work: the average 
ground-state energy of the bound state increases linearly 
with the concentration of impurities, however the slope is 
very small. Thus, for any reasonable concentration of im- 
purities we see that their effect on the bound state energy 
is very minimal. This is fully expected, given the analysis 
presented in the previous sections. Moreover, those re- 
sults also imply that increasing the relative strength a of 
the scattering potentials will not change this conclusion, 
since the effect of each scatterer saturates at large a. 

Finally, we comment briefly on the effects of other 
types of scatterers mentioned in the introduction. 
Clearly, attractive impurities cannot unbind the spin- 
polarized charge carrier - to the contrary, they will bind it 




impurity concentration (1/nm 2 ) 



FIG. 8: The ground state energy versus the impurity concen- 
tration. The energy was averaged over 1000 different impurity 
distributions and we chose a = 0.005. 



more strongly. Weak scatterers will have small, perturba- 
tional effects on the bound state. If the attractive poten- 
tial is extremely strong, a charge carrier becomes bound 
to the impurity itself and will screen its potential, mak- 
ing the scatterer "invisible" to the spin-polarized charge 
carrier trapped by the Zeeman potential. We conclude 
that such scatterers cannot damage the bound state. 

The other potential source of scattering is due to 
"noise" in the values of g e g reflecting the variations in 
the local density of Mn. One expects the typical scale for 
such variations (especially in digitally doped layers) to be 
very short compared to the tens of nm length scale of the 
bound wavefunction. If we model it in terms of a sum 
of short-range (delta function) noise, the conclusions will 
be as above, especially since the typical variation must 
be a small fraction of the average g c &, i.e. we expect this 
to be in a weak-scatterer regime. 



VI. CONCLUSIONS 

In this article we studied the effect of single and mul- 
tiple impurities on the spin-polarized state bound by a 
trapping potential generated by an inhomogeneous mag- 
netic field in a diluted magnetic semiconductor QW. We 
demonstrated that this effect is very limited, because of 
the large difference in length scales between the size of 
the wavefunction and the typical size of a strong scatter- 
ing potential. For very large scattering potentials, each 
impurity pokes a hole in the bound wavefunction, at a 
very small energetic cost. As a result, even for large con- 
centrations of very strong scatterers, the overall effect on 
the bound state is very limited. 

This result has obvious positive implications for any 
applications based on such trapped states, since they 
show that it is not necessary to worry about the effects 
of impurities or to use expensive methods to eliminate 
them. Samples grown in reasonably clean conditions 
should suffice for devices based on these spin-polarized, 
Zeeman-trapped states. 
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APPENDIX A: B-SPLINES 

B-splines have long been employed successfully in 
atomic and molecular physics^ but are not as freqently 
used in condensed matter physics^ The aim of this ap- 
pendix is to give a brief intuitive picture of the B-splines 
and their usefulness for solving Schrodinger equations 
with complicated potentials. Many more details on B- 
splines and their uses are available in the literature^ 

B-splines are piecewise polynomials and therefore are 
well suited for interpolation, having been extensively 
used in fitting tools, including many commercial software 
applications^ In the context of interest to us, they are 
used as a basis in which to expand the eigenfunctions. Of 
course, there are many possible finite basis sets and finite 
elements methods that can be used to solve a Schrodinger 
equations. The main advantage of using B-splincs is the 
flexibility to choose the grid points on which the B-splines 
arc defined. If we need to describe slowly-varying func- 
tions, a large mesh suffices, resulting in a small basis set. 
In regions where there are fast variations, one can use a 
finer local mesh, optimized to give the desired accuracy 
for the minimum increase in the number of basis func- 
tions. Furthermore, because the B-splines are piecewise 
polynomials, matrix elements can be efficiently evaluated 
to machine accuracy with Gaussian integration. Finally, 
the banded nature and sparsity of the resulting matrices 
allows for the use of very large basis sets, if need be. 

Suppose that we need to approximate a ID function 
in a given interval x £ [a, b]. We first define a knot- 
sequence of points in this interval {xi\a — xq < x\ < 
X2 < ... < xn =6}. The location of the points as well as 
their number will be chosen so as to optimize the process. 
For instance, the knot-sequence can consist of regions of 
equally spaced knots combined with regions where the 
knots are space in an exponential fashion, or any other 
suitable scheme. On this knot-sequence, we define the 
normalized B-splines of rank k by the following recursion 
relation: 

' w [ 0, otherwise v ' 



and 

Bi, k {x) = — - — ^ — S iife _i(a:)H — Xl+k ^ x B i+ i tk -i(x) 

where i is the index of a knot point and k is the order of 
the spline. Thus, the Bi^{x) is a polynomial of order k—1 
defined piece-wise in the interval [xi, Xi+k], an d vanishing 
outside it. Moreover, all derivatives up to the order of 
k — 2 are also continuous. Since we use these functions to 
expand eigenstates, which are continuous and have first 
and second continuous derivatives, it follows we should 
use a cubic (k = 4) or higher order B-splines. Here we 
use cubic splines, as they are the simplest splines with 
the desired properties. From now on, we simplify our 
notation and use Bi(x) to mean the cubic B-spline. The 
function of interest is then expanded in terms of these 
splines: 

f{x) = Y t a i B i {x) (A3) 

i 

where a are the N + 3 coefficients of our expansion (for 
the N + 1-point knot sequence). 

To solve an eigenvalue problem, we expand the wave- 
function in terms of B-splines and reduce the problem to 
a general matrix system: 

H eSv = 0, (A4) 

where H is the Hamiltonian matrix, S is the overlap ma- 
trix and v is the eigenvector containing the unknown co- 
efficients for the various B-splines. For a ID Schrodinger 
equation with a potential V(x), the matrix elements of 
H and S are equal to: 

^.^^^ 

Si,j = J Bi{x)B 3 {x)dx (A5) 

where we used an integration by parts in the kinetic en- 
ergy. The integrals for the kinetic energy and the overlap 
can be evaluated analytically, while those of the potential 
may need numerical integration, if V(x) is a complicated 
function. However, given the finite support of each B- 
spline, such matrix elements are non-vanishing only if 
\i — j\ < k, so the number of needed integrals scales like 
the number of basis functions, not like its square, as is 
the case for most other basis. 

The choice of the knot sequence plays an important 
role in these solutions. A good choice of knots distribu- 
tion can assure a fast convergence to the true eigenenergy 
with a small basis set. To illustrate this, we consider 
a gaussian potential and we calculate the ground-state 
eigenvalue and cigenfunction using two different knot se- 
quences. As illustrated in Fig [9], one is a uniform distri- 
bution of knots, while the other has an exponential dis- 
tribution of knots from the center of the potential. In the 
inset of Fig [SJ we compare the convergence of the energy 
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FIG. 9: The dots on the two lower lines illustrate the uni- 
form (squares) and non-uniform (diamonds) knot sequences. 
The curves show the estimated ground-state wavefunction for 
these N — 11 knot sequences, full line for the uniform and 
the dotted line for the non-uniform one. The inset shows the 
convergence of the energy to the true value as a function of 
the number N of knots used. 

as a function of the number of knots. It is clear that the 
non-uniform distribution is much more efficient in this 
case, giving a very accurate wavefunction and eigenen- 
ergy for a sequence with very few knots, i.e. a very small 
basis set. 

In the atomic physics, rigid boundaries are normally 
used when working with B-splines. Since in solid state 
physics we sometimes want to investigate bulk and trans- 
port properties, we opted to generalize the B-splinc 
method to periodic boundary conditions. The periodic 
conditions for B-spline are constructed in the following 
way (for cubic splines, although the generalization to 
higher orders is straightforward). Suppose we want B- 
spines on the support of the knot sequence {xi\a = xq < 

X\ < X2 < ■■■ < Xn = b}, 

i) Expand the knot sequence {xi\ to {x^} by adding 
the knots: 

x'_ 2 = £at-2 — (b — a) (A6) 
x_ x = XN-i — (b~a). (A7) 



Note that since x'_ x and x'_ 2 < a, they are outside the 
interval of interest. 

ii) Construct the B-splines {B®} with the usual proce- 
dure on the new support {x'j}. Note that the first two 
B-splines, which we call B^_^,B°_ 2 have finite support 
outside our desired interval. 

iii) Now construct the B-splines {Bi} by moving the 
pieces of defined outside the interval of interest, to 
inside [a, b] by a translation of period (b — a). 

Following the above procedure, all {Bi(x)} are periodic 
functions with period (b — a). Hence the functions in the 
Hilbert space expanded by {Bi(x)} with coefficients cjs: 



f(x) = J2 c i B i( x ) (A8) 

i 

are all periodic, thus ensuring the periodic boundary con- 
dition. Of course, for our problem of interest we choose 
b — a to be large enough that the localized wavefunction 
is all fully contained inside it. In other words, increasing 
b — a has no effect on the eigenenergies we calculate. 

The generalization to two-dimensional or higher- 
dimensional systems is straightforward. The wavefunc- 
tion needs to be expanded in products of B-splines 



ip(x,y) = ^a i , j B i (x)B j (y) (A9) 



and all the procedures above can be repeated. For our 
problem, we simply choose a high-density of knots in 
the immediate neighborhood of each impurity, where the 
wavefunctions change rapidly, and a wide, exponential 
mesh everywhere else, where the function changes slowly. 
The mesh near each scatterer has been optimized till con- 
vergence is reached. 
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